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Abstract 

This paper deals with spectral approximations for exterior elliptic 
problems in two dimensions. As in the conventional finite difference or 
finite element methods, it is found that the accuracy of the numerical 
solutions is limited by the order of the numerical farfield conditions. We 
introduce a spectral boundary treatment at infinity, which is compatible with 
the "infinite order" interior spectral scheme. Computational results are 
presented to demonstrate the spectral accuracy attainable. Although we deal 

with a simple Laplace problem throughout the paper, our analysis covers more 
complex and general cases. 
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1. INTRODUCTION 


In this paper we address some questions concerning applications 
of spectral methods to elliptic equations in exterior domains. Our 
goal is to investigate if one can obtain spectral accuracy in this class 
of problems. It is known through other methods that the accuracy of 
numerical solutions is governed by the order of the farfield conditions 
in addition to the order of the scheme itself. As an example for 
finite element methods see Bayliss et al. [BGT]. On the other hand 
boundary integral equations are very effective for exterior problems. 
They use appropriate Green’s functions which automatically take 
care of the farfield behavior (see [HMC] and [GW]). However, they 
are limited by the need to use explicitly known kernels and thus 
cannot treat variable coefficient problems. 

We believe spectral methods axe useful alternative ways to ad- 
dress these difficulties. Even more than in other numerical proce- 
dures, the formulation of the farfield boundary conditions is an es^ 
sential issue: a poor radiation condition may waste the high precision 
of spectral methods. A quite natural way of treating this problem 
within a spectral context is presented here, and forms the crucial 
part of this paper. 

The basic problem to be treated here is as follows: let D be a 
simply connected bounded domain in R 2 , whose smooth boundary 
will be denoted by T and whose exterior region by fi = R 2 - D . We 
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want to solve the problem 



Lu = 0 

in 0 

(1.1) 

u = g 

on r 


Boo(u) -* 0 

as r = 


where L is a second order uniformly elliptic differential operator in 
fi , g is a smooth data on T and the last condition represents the 
behavior of the solution at infinity (a radiation condition). 


As an example of a physical problem described by this formula- 
tion, one can think of D as a conductor in the field of a line source 
(located at * = Zq). Then L is the Laplacian operator A, u is the 
electrostatic potential, g(x ) = - log | * — 2o | an( ^ -®oo( u ) = u ~ 
log r. Another example is given by the incompressible, irrotational 
flow around a body D: again L = A, u is now the velocity potential, 
g(x) = 0 and Boo{u) = u - U 0 x - ^ log r, where U 0 is the main 
stream velocity in the x-direction and T is the circulation. Although 
(1.1) is the simplest model one can consider, we will see that the 
treatment by spectral methods is quite general. Even for eigenvalue 
problems - such as the Helmholtz equation, which is not solved here 
- the numerical farfield algorithm remains applicable. 


When the elliptic equation is discretized, the computational do- 
main is obtained by placing an artificial boundary Too surrounding 
the body at a finite distance. The radiation condition at infinity in 
(1.1) must be replaced by a boundary condition 


( 1 . 2 ) 


B(u) = 0 on Too, 
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which makes the problem mathematically well-posed and mimics the 
radiation condition at infinity. Such farfield conditions may be of 
local (or differential) type (see for instance [BGT], [G], [KM]), or 
of global (or integro-differential) type (see, e.g., [ADK], [FM], [M], 
[MCMj). They are usually obtained by approximately similar con- 
ditions satisfied exactly by the physical solution on Too . Hence 
they introduce an error on the numerical solution, which in principle 
should be comparable with the discretization error of the numer- 
ical scheme. The higher the order of the scheme, the higher the 
“order” of the farfield radiation condition. Since spectral methods 
produce discretization errors which decay faster than algebraically in 
the mesh-size, the farfield condition to be prescribed in conjunction 
with such methods must be particularly accurate to preserve the high 
precision of the interior scheme. 

The radiation condition we present here meets optimally this 
stringent requirement, because its error is just the truncation error. 
Our farfield condition is somewhat analogous to the global boundary 
condition used in [MCM], but it takes most advantage from being 
implemented in a spectral context, both in terms of computational 
efficiency and accuracy. Numerical evidence shows that our treat- 
ment of the radiation condition produces overall spectral accuracy 
on the solution. 

The plan of the paper is as follows. Section 2 is devoted to 
the discussion of the farfield conditions on the artificial boundary. In 
Section 3 we describe a pseudospectral algorithm for solving problem 
(l.l). Finally, Section 4 contains the results of some numerical tests. 
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2. AN INFINITE-ORDER RADIATION CONDITION 

Let us assume that the problem we want to solve is the electro- 
static potential problem, namely 

A u = 0 in fi 

(2.1) u = g on T 

u — log r bounded as r — * oo. 

Although we are dealing with the Laplacian operator, the reader 
should keep in mind that the spectral procedure we are going to 
describe is particularly designed for those situations - such as the 
case of operators with variable coefficients approaching a constant 
state at infinity - where the integral equation method is not feasible. 

The starting point for deriving several families of radiation con- 
ditions is the expansion of the exact solution into a convergent se- 
ries of eigenfunctions, usually through separation of variables in a 
neighborhood of infinity, where the differential operator has constant 
coefficients. For (2.1) we have 

(2.2) «(r, p) = log r + 

(r, tp) being the polar coordinates in the plane. Note that the right- 
hand side satisfies the radiation condition at oo. The coefficients 
a,/, are unknown. The farfield conditions are obtained by eliminating 
these constants, or a finite number of them, on the artificial boundary 


Bayliss, Gunzburger and Turkel use in [BGT] differential oper- 
ators in this process. Let us briefly recall their far-field conditions, 
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which we also implemented in a spectral context. The idea is to dif- 
ferentiate (2.2) m times in the r direction and eliminate the a*’ s for 
| k |< m through a linear cpmbination of such derivatives. For each 
m > 1 , this yields a differential operator 



which exactly anihilate the terms of order up to 2m in the series 
(2.2), i.e., 


(2.4) B m (u - log r - a 0 ) = 0 




, m = 1,2,, 


with ao estimated by averaging « in the farfield. The approximate 
solution u ap is required to satisfy the farfield condition 


(2.5) B m {u ap - log r - ao) = 0 on Too 

By comparing (2.5) with (2.4), it is seen that this method produces 
an ”a priori” (i.e., independent of the numerical discretization) error 
on the approximate solution, due to having dropped the terms on 
the right-hand side of (2.4). This error decays algebraically with the 
distance of the artificial boundary. It can be made arbitrarily small 
by raising the order of the farfield operator, but this may lead to a 
cumbersome or inefficient numerical procedure. 

The alternative approach which we follow in our construction of 
an infinitely accurate boundary operator consists of expressing each 
coefficient a ^ as a functional of u, rather than eliminating a finite 
number of them via a differential operator. Observe that for any 
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r > 0 , / rl fc l is the ft-th Fourier coefficient of the periodic function 

<p u(r, <p), as (2.2) shows. Thus the following representation holds 



If we differentiate (2.2) with respect to r: 


Mr,p) = ±[l-£|*|^ r e i ‘‘’ 

L k 

and use (2.6), we get an integro-differential relation satisfied by the 
exact solution of (2.1) on every circle of radius r. Precisely we have 

1 i 7r 

«r(r,p) = - 1 - — £ I k I e ik ^~^u[r,e)d6 

or 

(2.7) u T = -~- K*u, 

r r 

where * denotes the convolution of u with the singular kernel 

J oo 

K{q) = — y^kcoski]. 

We impose the radiation condition (2.7) on the approximate solution 
over the artificial boundary Too. Clearly, the precision of this farfield 
condition, as well as the efficiency of its implementation, rely upon 
the accuracy and the easiness with which the singular integral in (2.7) 
is evaluated. Accuracy and efficiency are guaranteed if the artificial 
boundary is a circle, which is not at all a restriction, and if the 
approximate solution is a trigonometric polynomial on Too, which is 
the case if a spectral Fourier method is used in the angular direction. 
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Assume that Too = {(r, ip) |,| r | = i?oo} and that the approx- 
imate solution, which we denote by u N is represented on by a 
trigonometric polynomial of degree N 

u n {R 0 o,<p)= bk [Roo)e tklp 

\k\<N 

(with ujv(i2oo) = u^y (f?oo))- Then one can easily check that 
K*u N = Yi \k\Z%(Rcc)e‘ kv , 

\k\<N 

i.e., the integral operator produces a new polynomial of degree iV, 
whose Fourier coefficients are obtained from those of w Ar (i? O o,0 sim- 
ply by multiplication by the modulus of the wavenumber. If u N is 
known on Too through its values u N (R co , <pj) at the equally spaced 
points jir/N,j = 0, 1, ...,21V — 1, rather than through its Fourier co- 
efficients, K * u N can be computed at the same modes exactly and 
efficiently, by Fourier transforming the values of u N to get its coef- 
ficients, then multiplying by the modulus of the wavenumbers and 
finally Fourier transforming back to get the point values of K * u N . 
Thus K * u N can be evaluated exactly in order lVlog 2 lV operations 
using a Fast Fourier transform, which is available anyway as part 
of the spectral calculations machinary. Moreover, the global char- 
acter of the convolution K * u N , which would destroy the locality 
of a finite difference or finite element discretization method, is not 
a drawback in a spectral context. Therefore, this kind of boundary 
operator appears to be very natural for our purposes. 

Thus the spectral solution is required to satisfy the radiation 
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condition 


( 2 . 8 ) 


it? = — 

r R 


OO 


1 — K * u 


N 


on Too. 


Unlike the family of farfield conditions proposed in [BGT], this is pre- 
cisely the same boundary condition satisfied by the exact solution, 
i.e., except for the truncation error (a finite number of modes), no er- 
ror is introduced on the numerical solution by the radiation operator 
on the artificial boundary. 


3. THE CHEB Y SHE V-FOURIER METHOD 


In this Section, we describe a pseudospectral method for solving 
the exterior problem (2.1), which uses the global radiation condition 
previously introduced. This algorithm consists of three main steps: 

i) mapping the physical domain onto the computational do- 
main; 

ii) building up the spectral collocation operator, which in- 
cludes the radiation condition; 

iii) solving the resulting system of discrete equations by var- 
ious relaxation techniques. 

Let us describe each step separately. 

i) The Coordinate Transformation 

It is not restrictive to put the origin of the coordinates within 
the domain D. We assume that the boundary T of D is represented 
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in polar coordinates by the equation r = R{<p), 0 < <p < 2 ir, where R 
is a smooth function. The exterior problem (2.1) will be discretized 
in the anular region (l between F and the artificial boundary F©©, 
which is a circle of radius JZ©© : 

n = {(rcosp, rsin ip) eR 2 | 0 < <p < 2 jr, R ( <p ) < r < JZ©©}. 

This domain is mapped smoothly onto the computational domain, 
the rectangle 

n = {(s, 6) e R 2 ) ) — 1 < s < 1, 0 < 0 < 2ir} 

in such a way that the image of T is the side s = — 1. Let us remark 
at the outset that we do not attempt to obtain any special proper- 
ties in the map (r, <p) — * («, 0), except, of course, boundary fitting. 
In particular, the map is not conformal, as this would produce an 
immediate solution. 


The coordinate change has the form 


(3.1) 


ip = 0 

r = r(«, 0) 


with r(-l, 0) = iZ(0) and r(l, 0) = iZ©©. 0 < 0 < 2ir. The func- 
tion r(«, 0) may be chosen to be an exponential stretching in the 
s-direction, precisely 


(3.2) r(s,0) = iZ(0) exp (a(0)(s + l)) 

with a(0) = ^ log (iZoo//Z(0)). At least two reasons suggest this 
choice. First, the exact solution of (2.1) satisfies for r large enough 


u log r = log(iZ(0)e a ^ ff+1 )) c* a(0)s. 
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Accordingly, u at infinity is a power of s, which is correctly approxi- 
mated by Chebyshev polynomials in s. The next reason is the expo- 
nential stretching avoids the unnecessary crowding near the artificial 
boundary of the Chebyshev points involved in the colloation process 
(see [G] for similar considerations). 

Numerically, the exponential stretching produces more accurate 
results than other stretching functions of polynomial type (see [B])- 
However, for some classes of problems algebraic stretching is the best 
(see [B], [GO]). The choice is dependent upon the problem. 

ii) The Pseudospectral Operator . 

When written in the (s, 0) coordinates, the Laplace equation 
takes the form 

(3*3) Lu = au ss + bu s $ + cugg + du s = 0, 

for some variable coefficients a,b,c,d satisfying the ellipticity condi- 
tions a > 0, c > 0 and ac — b 2 > 0 in fi. The solution u is periodic in 
6, satisfies a Dirichlet condition at s = — 1 and the radiation condition 
(2.7) at 8 = 1. Therefore we look for an approximate solution u N 
which is a trigonometric polynomial of degree N in 0 and an algebraic 
polynomial of degree M in $: 

M 

«"(*.«) =E E T m (.) e M - 

m=o |/c|<jV 

In the sum the asterisk means that u% N = u*(-n) for a11 m > and 
T m (s) is the m-th Chebyshev polynomial of the first kind. u N 
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is equivalently determined by its values ujy at the (M + l)x2 N 
Chebyshev-Fourier nodes 

(*(,«>) = (cos 0, £), i = 0, .... M j = 0, .... 2JV - 1. 

While for j=0,...,21V-l are prescribed by the Dirichlet condition 
on r, the unknowns u fj for i=0,...,M-l, j=0,...,2iV-l are computed 
by solving the set of collocation equations 

(3.4) £,,(«*)(*, fy)=0 

The operator £ sp is pseudospectral approximating to the operator L 
defined in (3.3), which makes use of the radiation condition (2.8) - 
as well as of the Dirichlet condition on T - in the evaluation of the 
partial derivatives of u N at the mesh points. Let us describe this 
process. The derivative does not involve boundary conditions, 
hence it is computed by a standard procedure through the discrete 
Fourier transform. The derivative is computed spectrally from 
all the values of u N on the mesh. On the artificial boundary, what 
is obtained does not necessarily satisfy the radiation condition (2.8), 
which in the (s, 6) coordinates takes the form 

1 - K * u N 

Thus we modify the values of u ^ on the artificial boundary according 
to (3.5). The term K*u N is evaluated in the way described in Section 
2. We end up with a modified partial derivative of u N , let us call it 
which we further differentiate spectrally to produce u^ s and u^ $ . 

In this procedure, the boundary conditions are imposed implic- 
itly, namely during the evaluation of the spectral operator £ sp . Thus 


(3.5) 


N _ 

Us ~ dr 
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lisp is an affine operator, since the boundary conditions are non- 
homogeneous. In view of the further discussion, it is convenient to 
write (3.4) in the equivalent form 


(3.6) 


(L sp u N - /)(«,-, $ s ) =0 * = 0, ...,M - 1, 

j = 0,...,2JV-1 


where L sp denotes the linear part of £ 9p and the vector / takes into 
account the inhomogeneity in both the Dirichlet and the radiation 
boundary conditions. 


We conclude this section with a few remarks on the implemen- 
tation of the local boundary conditions (2.5). For m=l we have a 
Robin-type condition, which can be implemented implicitly as de- 
scribed previously, by modifying the computed values of on the 
artificial boundary to satisfy the radiation condition. For m>2,B m 
involves higher order radial operators. They may be computed spec- 
trally, and (2.5) can be imposed explicitly, as a set of linear equations 
at the points on Too (in this case the elliptic equation is collocated 
at the interior points only). However, the resulting algebraic system 
exhibits wild eigenvalues which does not compromise the possibility 
of using efficient iterative methods of solution. So, it is preferable to 
eliminate the higher order radial operators by means of the differen- 
tial equation (3.3) on the boundary. This process, although rather 
complex, produces a boundary operator which is just first order in the 
radial direction; the associated boundary conditions can be imposed 
implicitly, as described above. 
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iii) Iterative Methods for the Discrete Problem . 

The direct inversion of the algebraic system (3.6) is extremely 
impractical since the associated matrix, which we still call L sp , is full 
and severely ill-conditioned. This matrix is not even computed, as it 
would require too much storage. Instead, the matrix-vector multipli- 
cation L gp v can be done efficiently through fast Fourier transforms, a 
process which needs just few storage matrices of size Mx2N for u and 
its partial derivatives and the coefficients of the operator. Therefore 
we always resort to iterative methods of solution, and specifically 
to those methods which update the solution by few matrix-vector 
multiplications of the type L sp v. 

Explicit time advancing techniques are a source of such methods, 
as iteration may be thought of as an evolution towards steady state of 
the solution of ut = L sp u-f. Convergence requires the stability of the 
scheme for the operator of the problem; moreover, the convergence 
history is heavily influenced by the conditioning properties of the 
operator and by the choice of the pseudo-time step. 

Among these methods, we tested the DuFort-Frankel algorithm 
([G Li], [G L 2 ], [F]), in the form 

where a and c are the coefficients in (3.3) and Asy = sy_ i — Sj, Ad = 
7 r/jV. This algorithm is known to be unconditionally stable for con- 
stant a, c, As and for Dirichlet or Neumann boundary conditions, 
once <r has been taken large enough. However, for variable coef- 


a ij 






ok/ 1 -2< J+ t <?-■),} 


(3.7) 

n+1 


{ 




— U 


2A t 
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ficients and radiation boundary conditions there are limitations on 
At, and the scheme (3.7) was found to converge very slowly to steady 
state. Even when the parameters At, <r were fixed at optimal values 
(by trial and error, since there are no formulas for the general case), 
or dynamically redefined at each time step the convergence was still 
slow. 

Due to the ill-conditioning of spectral operators, it is generally 
accepted that iterative methods for the solution of steady state prob- 
lems have to be used together with preconditioning techniques. As 
first pointed out in [D] and [Mo], low order finite difference approx- 
imations to the boundary value problem at the collocation nodes 
provide a good way of preconditioning spectral methods. The matri- 
ces obtained by this process, denoted generically by Ljd are sparse 
and hence easily invertible. Moreover, the condition number of the 
matrix L L sp is bounded independently of the mesh size, or it grows 
slowly with it. These facts have been proved for some constant coef- 
ficient operators ([0],[HLAD]), and observed numerically for a wide 
class of variable coefficient operators, subject to Dirichlet or Neu- 
mann boundary conditions (see, e.g., [CQ],[HLAD],[PZH]). 

Hereafter, we will describe a finite difference preconditioning 
matrix Ljd for the spectral operator with integral radiation condition 
introduced in this section. At the mesh point (i,j) the operator (3.3) 
is approximated by second order finite differences over the molecule 
{(tj), {i ± 1 ,j), ( i,j ± 1), (i + 1 ,j + 1), (t — l,j — 1)}. This yields 
a 7-diagonal matrix. The points (i + l,j + 1) and (*' — 1, jr — 1) 
are needed for consistency with the mixed term u s q. A 5-diagonal 
matrix is obtained instead by approximating the reduced operator 
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a u S9 + cuee over the molecule {(»,/), (* ± l,j), (i , j ± 1)} The latter 
matrix can be inverted less expensively than the former one, and it 
provides a fairly good preconditioning as long as the body shape does 
not deviate too much from a circle. Otherwise, the term u s $ in the 
differential operator is significant, and the 7-diagonal matrix has to 
be preferred. 

In both cases, the outer points of the computational molecules 
centered on the artificial boundary are eliminated using the homo- 
geneous form of the radiation condition (2.8). The exact evaluation 
of the integral term K * u would produce a full 2Nx2N block cor- 
responding to the unknowns on the artificial boundary. Instead, the 
band structure is preserved if the global operator is approximated by 
a three-term formula, 

(K * u) 03 - c* k + u 0 ,j+i + k 0 u O j + k-Uoj-i. 

The coefficients k± and ko are determined in such a way that the 
approximation is exact for u = 1, cos 6 and sin 6 on the boundary. 

At the points on the body surface, homogeneous Dirichlet con- 
ditions are implemented. 

In this way we end up with a matrix W, which has 7 non- zero 
diagonals (5 if the five-point molecule is used), plus 2(M-1) non zero 
entries, due to the periodicity condition in the 6 direction. An in- 
complete LU decomposition of this matrix can be built up following 
Meijerink and van der Vorst [MV]. Precisely, a lower and an upper 
triangle matrix L and U are introduced, respectively with the same 
(a priori) non-zero elements as the lower and the upper triangular 
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parts of X. U is normalized to be identically 1 on the main diagonal. 
These matrices are defined by the condition that the product Mi 
equals X at the (a priori) non-zero elements in X (A different condi- 
tion has also been implemented, essentially with equivalent results: 
Mi and X must agree at the off-diagonal non- zero elements of X, 
and the row sums of Mi and X must be equal, see [W]). Finally, we 
set Lfd = Mi. This matrix is a close approximation to X, while it 
allows the inversion of a linear system in O(MN) operations. 

The algebraic system (3.6) with this preconditioning has been 
solved by Richardson iterations 

(3.8) u n+1 =u n +aLjt{f-L sp u n ) , n > 0. 

Alternatively, one can use Dufort-Frankel iterations 

(3.9) u n+1 = u n + - L sp u N ^j -i(u n+1 -2u n +u n-1 ), 

which corresponds to a second order method. The following discus- 
sion concerns (3.8), but the conclusions hold true for (3.9), too (see 

|CQ|). 

The choice of the acceleration parameter a in (3.8) is quite deli- 
cate. Formulas for a involving the eigenvalues of A = Lj^L sp (see for 
instance [O]) are impractical, since no theoretical information on the 
spectrum of A is known. Attempts to estimate dynamically the ex- 
treme eigenvalues of A (see [HLAD] for such an algorithm) have given 
poor results in our situation. Instead, a minimal residual strategy 
for choosing a at each iteration proves to be quite efficacious. Pre- 
viously, a n is determined in order to minimize the £ 2 -norm of the 


17 


residual r n+1 =f- L sp u n+1 , (see [WH]), i.e., 


(3.10) 


(• L 8p P n ,L ap p n Y 


where 


(3.11) p n =L J lr n . 

With this choice, (3.8) is nothing but a descent method of the form 

(3.12) u n+1 =u n +a n p n , 

where the descent direction p n is determined according to (3.11). 
This method is quite appealing, since it produces fast convergence 
with modest work and storage requirements. Its drawback lies in 
the fact that it may break down if the symmetric part of the ma- 
trix LgpLJt is indefinite, as a n may become zero. This occurrence 
has never been observed in the case of full Dirichlet boundary con- 
ditions, or constant coefficient elliptic operators. Unfortunately, the 
breakdown does happen for the problem under consideration, in cases 
where the geometry is far from being circular in shape. Then, a 
descent direction other than (3.11) must be used to continue the 
method. 

We found it very effective to determine p n according to the 
method Orthodir, proposed by Young and Jea ([YJ]). In its origi- 
nal formulation p n is computed to be A T A - orthogonal to all the 
previous descent directions (recall that A = Lj^L sp )^ ? 

71 — 1 

p n = Ap n ~ l - ^2 Pn*P k ’ 

k=0 


(3.13) 
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with 

« ■ _ UV- 1 , V s ) 

{Ap*,Apk) ' 

The method is guaranteed to converge even if the symmetric 
part of A is indefinite. We use a truncated version of the method, in 
which the sum in (3.13) starts from k=n - 2 (for this version, however, 
convergence is not assured anymore). One step of Orthodir requires 
twice as many operations as one step of the original method (3.8). 
Thus, we shift to Orthodir whenever the a n given by (3.10)-(3.11) 
becomes too small, say or” < .001. After one step of Orthodir, we 

continue with Richardson’s iterations (3.8). No break-down of this 
algorithm was observed. 
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4. NUMERICAL RESULTS . 

Evidence about the good performances of the method proposed 
here can come - at the moment - from numerical tests only. Indeed, 
the present status of the analysis of spectral methods is unable to 
handle such ” complex” situations as multidimensional, variable co- 
efficient equations with integrodifferential boundary conditions (see 
[GHV] for an overview of mathematical results about spectral meth- 
ods). In particular, analysis of variable-coefficient operators are 
known for boundary value problems of Dirichlet type only (see the 
paper by Canuto and Quarteroni in [GHV]). Results for the implicit 
imposition of boundary conditions of Neumann or third type have 
not been established yet, except for the one-dimensional, constant 
coefficient case. 

The hard point is to prove that the collocation scheme described 
in Section 3 is stable; namely, that the Chebyshev-Fourier solution 
to the problem 

AtV* = F, in a 

(4 1) = G ' onT 

+ i k,w"-H on Too. 
dr r 

satisfies an estimate of the form 

(4.2) || W N ||~< || F || ft + || G || r + || H ||roo }, 

(with C > 0 independent of N), where the norms in (4.2) axe norms 
of suitable Sobolev spaces. If we assume that (4.2) holds, then it is 
immediate to prove not only that the approximation is convergent, 
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but also that the error decays spectrally. Actually, if u is the solution 
of (2.1), set u = I M P N «, where P/v denotes the truncation of a 
trigonometric series in 6 to the order N and I M denotes the algebraic 
interpolation of degree M in the 5-direction at the collocation points. 
It is well known that the error between u and u decays spectrally. 
Since Pn[K *u) = K * P^u, it is readily seen that W N — u N — u is 
the collocation solution of the problem 


AW" = A(u - u) 

in n, 



W N =g-P N g 

on r 



dW N l 

—Z h -K * W N 

or r 

II 

Q?| Cd 
| e 

, du 

dr 

on 


All the right-hand sides in ( 4 . 3 ) decay faster than algebraically in the 
Sobolev norms as N, M — oo, hence so does W N , according to (4.2). 
We conclude that u N approximates u with spectral accuracy. Note 
that if ( 4 . 2 ) holds, the convergence is guaranteed for any fixed value 
Poo °f the radius of the artificial boundary. 

Let us present now some numerical results. 

In order to compare the behavior of our global radiation condi- 
tion, the Bayliss-Gunzburger-Turkel radiation conditions of the first 
and second order were also implemented. The same pseudospectral 
algorithm described previously was used. 

In order to appraise our numerical results we need exact solution 
of the problem for an arbitrarily shaped domain D. We simulated this 
situation by considering a point source at = {x 0 ,y 0 )e D, whose 
exact solution is 


log|z-z|, 
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We prescribed this value on T. Note that this satisfies (2.1) with 
oo = 0. The numerical procedure begins only with this value on T. 
Moreoever, one can add to this above solution any harmonic function 
which decays at infinity and still satisfied by (2.1). Such solutions 
are used for comparison purposes. We considered three significant 
situations. 


CASE 1: the body is a circle of radius 1, and the exact solution 


is u(x,y) = log 


(x - .7) 2 + y 2 




CASE 2: the body surface is described by the function R(^) = 

1 + .4 cos 2<p (see Fig. 1), and the exact solution is 

. , , cos 2 <p ^cos3op 

u(r, <p) = log r + — — + 3—t—. 


CASE 3: the body surface is described by the function R(<p ) = 

1 + .7 cos ip + .3 sin 2<p (see Fig. 2), while the exact solution is 


u(x,y) = log (x - l) 2 + [y- l) 2 


The algorithms were tested on a 9x8, a 17x16 and a 33x32 grid, and 
for different positions of the artificial boundary. If we define the rela- 
tive distance of the artificial boundary to be p = Roo / 
then the range of p was the interval [1.2,5.]. 


Table 1 provides some information about the behavior of the 
iterative algorithm used to solve the system (3.6). We report the 
relative residual (i.e., the l 2 - norm of the current residual divided 
by the l 2 - norm of the initial residual) at some selected iterations. 
The convergence histories are relative to the intermediate case p — 3, 
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and to the use of the integral radiation condition (2.8). These results 
can be considered typical, in the sense that the rate of convergence 
becomes just slightly worse where p approaches 1 (because of the in- 
creased stretching of the coordinates), and it is essentially the same 
if the local radiation conditions are imposed instead. It is observed 
that the preconditioning is sensitive to the geometry (i.e., to the 
coefficients of the differential operator). Moreover, the spectral ra- 
dius of the iteration matrix depends on the mesh-size (although some 
strange behaviors appear, such as in Case 2 where the convergence 
on the 17x16 grid is faster than on the 9x8 grid, or in Case 3 where 
the histories for N=16 and N=32 are essentially similar). These facts 
are clearly related to the use of an incomplete LU decomposition in 
the inversion of the finite difference system. However, the precon- 
ditioning is globally effective. The reader should take into account 
that over 1000 iterations of the Du-Fort-Frankel algorithm (3.7) are 
needed to drive the residual below 10“ 10 , in the simplest geometry 
of the circle. Moreover, in all the cases tested, at most some 50 it- 
erations of our algorithm were enough to produce 2 significant digits 
in the spectral solution. 

Tables 2-4 contain the results of our tests, which allow to com- 
pare the performance of the first order (F), the second order (S) 
[BGT] radiation conditions, and the integral (I) radiation condition 
on the issue of the spectral accuracy. In each cell, the upper number 
is the relative error between the exact and the spectral solution in the 
/ 2 - norm at the grid in the computational domain. This error may 
appear geometry-dependent, since the grid points are displaced in 
the physical domain by the stretching of coordinates. For this reason 


23 


we have included the relative error in the maximum norm, which is 
independent of the mapping of the domain (this is the lower number 
in the cell). The two errors behave qualitatively in the same way. 

In general, it is seen that the integral radiation condition always 
produces the best results, and it is the sole to guarantee the spectral 
accuracy in all the cases. 

This is evident in Case 1 (Table 2), where the I-condition pro- 
duces spectral accuracy for any value of p. Thus the error just comes 
from the discretization scheme, while it is weakly sensitive to the po- 
sition of Too. On the contrary, the F- and the S- conditions produce 
spectral accuracy only when Too is far enough from the body. Other- 
wise, the truncation error introduced by the local radiation condition 
dominates the discretization error, and the global error increases as 
Too is brought close to the body. This is particularly evident on the 
33x32 grid. The reason is that in the expansion (2.2) of the exact 
solution, the terms with index | k | larger than 2 (or 3) are treated 
incorrectly by the local radiation conditons, and correctly by our 
global condition. 

This phenomenon is even clearer in Case 2, where the exact 
solution is chosen on purpose to exhibit the different behavior of the 
radiation conditions. The I-condition is the only one which allows the 
artificial boundary to be brought very close to the body, yielding the 
smallest errors. Note that on the 9x8 and the 17x16 grids, the error 
decreases as Too approaches D even when the F- or the S- conditions 
are used: the error in this case is dominated by the discretization 
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error, which takes advantage by the refinement of the mesh. But on 
the 33x32 grid, the same behavior as in Case 1 occurs. 

Finally, Case 3 is an example of a situation where an exact solu- 
tion with a singularity relatively close to the boundary of the body is 
coupled with a fairly complex geometry. Again, the I-condition per- 
forms the best, producing spectral accuracy. It should be noticed, 
however, that now even the error obtained by the I-condition gets 
worse if the artificial boundary is too close to the body. Since the 
exact solution is qualitatively comparable with the exact solution of 
Case 1, the origin of this phenomenon should be found in the fact that 
the coordinate mapping produces larger variations in the coefficients 
of the elliptic operator, hence there is an overall loss of precision. 
This remark may suggest the use of a multidomain technique, when- 
ever the geometry is too complex to be handled accurately by the 
one domain mapping discussed here. 
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Figure 1. 
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Figure 2. 
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TABLE 1: Convergence histories (* means: below 10 ) 
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